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We propose a classical SU(2) gauge field in a fiavor-space locked configuration as a species of 
radiation in the early Universe, and show that it would have a significant imprint on a primordial 
stochastic gravitational wave spectrum. In the flavor-space locked configuration, the electric and 
magnetic fields of each flavor are parallel and mutually orthogonal to other flavors, with isotropic and 
homogeneous stress energy. Due to the non-Abelian coupling, the gauge field breaks the symmetry 
between left- and right-circularly polarized gravitational waves. This broken chiral symmetry results 
in a unique signal: nonzero cross-correlation of the cosmic microwave background temperature and 
polarization, TB and EB, both of which should be zero in the standard, chiral symmetric case. We 
forecast the ability of current and future cosmic microwave background experiments to constrain 
this model. Furthermore, a wide range of behavior is shown to emerge, depending on the gauge field 
coupling, abundance, and allocation into electric and magnetic field energy density. The fluctuation 
power of primordial gravitational waves oscillates back and forth into fluctuations of the gauge field. 

In certain cases, the gravitational wave spectrum is shown to be suppressed or amplified by up to 
an order of magnitude depending on the initial conditions of the gauge field. 
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The recent flourish of attention to the imprint of grav¬ 
itational waves on the cosmic microwave background 
(GMB) [1, 2] has spawned an array of models that af¬ 
fect the physics of these tensor modes. On the obser¬ 
vational side, the leasing R-mode spectrum of the GMB 
has been measured to ever increasing accuracy and B 
mode surveys are getting close to detecting the primor¬ 
dial spectrum [2-8]. This polarization information will 
put stringent constraints on theories of the early Uni¬ 
verse. However, the subsequent propagation of tensor 
modes is also affected by the composition of the Universe 
which leaves an imprint on the tensor spectrum [9-12]. 
The prospect of precision measurements of cosmological 
tensor modes opens up this area of research extending 
the era of precision cosmology to more intricate models 
concerning late-time tensor propagation. 

A particular class of ideas proposes cosmological parity 
rotation (GPR). The standard model of particle physics 
is invariant under simultaneous reversals of charge (G), 
parity (P) and time (T). However, it does not obey par¬ 
ity invariance on its own, most famously demonstrated 
for weak interactions [13]. This motivates similar be¬ 
havior in cosmological physics and CPR is an example 
thereof. Among other symptoms CPR affects proper¬ 
ties of gravitational waves as they possess intrinsic par¬ 
ity. First and foremost, allowing for different behavior 
for left- and right handed tensor modes results in GMB 
cross-correlations that are typically not allowed. Some 
mechanisms that cause this physics have their roots in 
novel interactions of electromagnetism, as for example 
cosmological birefringence [14], which rotates E into B- 
modes. On the other hand the gravity sector itself can 
source the parity violation through explicit changes in 
the action. A wide variety of theories has been devel¬ 
oped that accommodates these effects [15-17]. 

Parity violation happens naturally in the weak sector 
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of the standard model of particle physics. However, the 
standard model is presumably only a low-energy limit 
of a grand unified theory. In cosmology both inflation 
and late-time cosmic acceleration require new physics be¬ 
yond general relativity - and from a particle physics point 
of view - beyond vanilla ACDM. It is therefore natural 
to ask whether the new physics responsible for cosmic 
inflation or late-time acceleration incorporates parity- 
violating effects. 

Inflation is usually thought of being caused by a scalar 
field. However, there has been a growing number of vec¬ 
tor and gauge field models appearing in the literature 
over the past decades [18-29]. In this paper we take in¬ 
spiration from vector field inflation and Yang-Mills (YM) 
theories to introduce a fluid that causes parity-violating 
effects in the Universe. Field configurations that render 
the YM field homogeneous and isotropic have been found 
and these are employed to hold up isotropy in the Uni¬ 
verse. We propose adding a dark YM radiation fluid to 
the energy budget of the Universe. This YM fluid violates 
parity on cosmological scales and changes the behavior 
of gravitational waves significantly. We use this as a sim¬ 
ple model for chiral effects in gravity without requiring 
changing the Einstein-Hilbert action. 

In this paper we will show the effect of the YM fluid 
on cosmological observables. We compute the change 
in gravitational wave amplitudes and their correspond¬ 
ing energy densities. We determine how this causes de¬ 
viations from vanilla ACDM in the CMB spectra and 
conclude by computing forecasted constraints that mea¬ 
surements of the classically vanishing cross-correlations 
impose. 

The paper is structured as follows: We introduce the 
model in Sec. H and derive the perturbation equations 
for the scalar, vector and tensor modes in Sec. HI. We 
solve these equations in Sec. IV where we separate off the 
treatment of the CMB into Sec. V. In Sec. VI we derive 
constraints that experiments put on this model and we 
conclude in Sec. VH. Details about the implementation of 
this model in CAMB [30] and solutions to the background 
equation of motion can be found in the appendices. 


II. BACKGROUND MODEL 

We consider the standard cosmological model with the 
sole addition of a new gauge field as a toy model for 
cosmological chirality violation. This effect is not equiv¬ 
alent to parity violation in fundamental physics where 
some terms in the Lagrangian violate parity explicitly. 
In spirit, the introduction of a particular configuration of 
the gauge field rather acts like the vacuum expectation 
value of a scalar field spontaneously breaking symmetry. 
Here, the broken symmetry is chirality: Left- and right 
handed gravitational waves will behave differently upon 
introducing this field. Additionally, we treat the gauge 
field as a classical, low-energy limit of a YM gauge field. 


The action of the theory is given by 

S = J + (1) 

where Greek letters are used to represent space-time in¬ 
dices, lowercase Latin letters are spatial indices, and up¬ 
percase Latin letters I are reserved for the gauge group 
indices. External matter is included in the Lagrangian 
density Lm- The constant Mp^ = SttG is the reduced 
Planck mass. The kinetic term for the gauge field is 
canonical and clearly inspired by high-energy physics, 
specifically Yang-Mills gauge theory. Symmetry break¬ 
ing is a common feature of these theories and we use 
similar interactions to model chiral effects in cosmology. 
Therefore the field strength tensor is taken directly 
from Yang-Mills-like theories and it reads 

Ffiv = dfiAl - Ajfj^AKiy ( 2 ) 

where ^ym denotes the group coupling constant. In the 
internal space indices get raised and lowered with the flat 
metric Sj. 

The fluid equations for the gauge field follow from vary¬ 
ing the action (1) with respect to A^: 

= 9^"kA>^jF^, (3) 

and the stress energy tensor for the gauge field is 

(4) 

Many authors have used gauge fields in cosmology - for 
inflation as well as for late-time behavior [18-29]. A com¬ 
mon issue is that spin-1 fields have a preferred direction 
and therefore break isotropy. Nevertheless, there is a 
mechanism that provides isotropic stress energy [22, 31] 
which relies on the homomorphism between the group 
of space rotations, 0(3), and the internal group space 
SU(2). Making this choice for the gauge group, too, al¬ 
lows for isotropy: aligning the global part of the gauge 
SU(2) with the rotational SU(2) restores rotation sym¬ 
metry because gauge fields are only defined up to these 
gauge transformations [32]. Therefore I S {1,2,3} for 
each of the SU(2) generators. 

We consider cosmological solutions wherein directions 
of the internal SU(2) space are aligned with the princi¬ 
ple axes of the Cartesian, spatially flat Robertson-Walker 
space-time, = —dt'^ + a^[t)dx^. Specifically, we as¬ 
sume a flavor-space locked configuration for the gauge 
field, wherein 

imsf f = i ,5, 

[0 g = 0 ^ ’ 

with all other components vanishing where the held (j){t) 
is a homogenous scalar. This particular conhguration 
is just an example for the group of assignments of fit) 
to the gauge held that is spanned by all SU(2) gauge 
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transformations on Aj^. In particular one can rotate the 
^(t) scalar inside without affecting the background 
or perturbation equations. This just follows from gauge- 
invariance of the YM Lagrangian under SU(2) transfor¬ 
mations Aj^ ^ Aj^ + where 0^ is a 

local, infinitesimal SU(2) gauge transformation. For the 
rnonzeroest of this analysis we choose the mapping in 
Eqn.(5). 

With this assignment the nonzero components of the 
field strength tensor are 

^oi = and Flj = ij 

where a dot denotes a derivatives with respect to cosmic 
time. This field configuration resembles a condensate of 
massless gauge bosons, like the electroweak bosons before 
symmetry breaking, with their three hypercharge spin 
vectors pointing along the x, y and z directions. The 
Lagrangian expressed in terms of the scalar cj) is 

= . ( 6 ) 

Using the above equations, from Eqn. (4) the energy den¬ 
sity and pressure of the gauge field are 




I 2 

+ gvM 





I 2 

+ gvM 



(7) 


therefore, the standard, massless SU(2) gauge field has 
an equation of state of w = 1/3, like radiation. The 
energy density in Eqn. (7) is composed of ‘electric’ (j) and 
‘magnetic’ contributions, motivated by the term 

structure of the gauge field. 

To gain some insight into the dynamics, we isolate the 
YM sector of the model. We write the gauge field action 
in the spatially flat Robertson-Walker coordinates and 
plug in the above expressions, Eqn. (6), but then switch 
to conformal time whereby 

S = J dtd^x^ 

= J dr d^x^ - gyM^fj)'^) , 

and derivatives with respect to conformal time are de¬ 
noted with a prime. The equation of motion for this field 
in conformal time is 

4>" + 2gYM^/>^ = 0. (8) 

Interestingly, this description in terms of conformal time 
is independent of the scale factor, which allows for a 
general, analytic exploration of the properties of this 
field. Some details are given later in this section and 
Appendix B. 


The solution may be expressed in the form of a Jacobi 
elliptic function [33] 

5YM/>('r) = Cl sn(ci(T - n) + C 2 \ - I). (9) 


The constants are determined at the initial time as 

C2 = F(csc-i (I + I - I), 

where F is an elliptic integral of the first kind. As shown 
in Appendix B, (/ is an oscillating function of time. The 
smaller the value of ^ym , the slower the rate of oscillation. 
However, for any coupling, the field behaves as radiation 
with equation of state w = 1/3. 

We rephrase the constants Ci and C 2 in terms of a pa¬ 
rameter describing the energy density in the field, and a 
parameter describing the allocation of energy in the elec¬ 
tric and magnetic components. The energy density in the 
gauge field as a fraction of the photon energy density is 
RyM = Pyu/P-y so that the initial values of the field and 
its derivative are 


(j)\ — FIqMp^ 2i?YM^7,o sin^ 9 

,2 HoMp -2. 

(j>i = - \ zAymS 0 0 cos^ 0 

gym ^ 

where 9 S [0, |] dials between electric and magnetic field 
energy. It will be useful to introduce two inverse comov¬ 
ing length scales 


fco = fiRyyiHoaQ 

kg = \J gyuMpaoko 


( 12 ) 


so that Cl = kg and C 2 = F(csc ^(•\/sec0)| — 1) in 
Eqn. (10). The gauge field oscillates with period 


V^kg ’ 

but its energy density and pressure scale with equation 
of state w = 1/3, like radiation. For the field to remain 
coherent on cosmological time scales (and avoid a secular 
instability in the linear perturbations), and thereby have 
the maximum effect on cosmological physics, the cou¬ 
pling must be exponentially small, ^ym ~ 0{Hq/Mp) ^ 
10“®°. This small number can be achieved if our theory 
originates with a dilatonlike factor e°'F^. With this cou¬ 
pling a dilaton can roll high up without additional cost 
since vanishes under equipartition of energy between 
electric and magnetic modes. Supposing that a mecha¬ 
nism stabilizes a there would be no effect on the equa¬ 
tions of motion, but the gravitation of the gauge field 
stress energy tensor would be magnified by this factor. 
Since our requirement that the field remains coherent on 
cosmological time scales is equivalent to setting kg < Hq, 
then an original coupling g ~ 0(1) can be engineered if a 
suitably large value of a is permitted. An origin for this 
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small coupling as well as the flavor-space locked config¬ 
uration might be devised in an inflationary epoch along 
the lines of [34, 35], although that is beyond the scope 
of our present investigation. Inflationary scenarios based 
on a similar gauge field that include self-interactions with 
couplings to matter fields, have been studied elsewhere 
[36, 37], 

Other potential constraints on this model could come 
from electric dipole measurements (EDM) assuming that 
baryons and leptons are charged under the gauge group of 
this field. Generally, EDMs violate both parity and time 
symmetries and they yield model-independent measures 
of CP violation in nature. However, as we will show now, 
the effects of the YM field on such a dipole relaxation pro¬ 
cess are negligible. We assume that a particle of mass M 
possesses an intrinsic dipole moment ^ = ( 7 ym/ 2 M under 
this gauge group. The spin-flip energy due to the cou¬ 
pling between the dipole moment and the flavor-space 
locked field is therefore 2/iH where B = gyu'p'^- Us¬ 
ing B - MpHo^RYu/Zeq where Zgq is the redshift of 
radiation-matter equality, and 5 'ym ~ Hq/Mp, then the 
energy shift is roughly {Hq/M)^R/Z eq. Using these re¬ 
sults, we see that particles with mass M produce photons 
with a wavelength ~ Hq^{M/Ho), which for any reason¬ 
able mass is many orders of magnitude longer than the 
Hubble horizon scale today. 


III. PERTURBATION EQUATIONS 


In this section we study the most general perturbations 
of the background solution presented in Sec. H. We are 
interested in linear perturbations in this paper. These 
can be split up into scalar, vector and tensor perturba¬ 
tions that are decoupled from each other. 

We are following Subsection HI A of [38] for develop¬ 
ing the perturbation quantities. In general there are 9 
perturbation degrees of freedom for the gauge field and 
6 in the metric sector. Some of the perturbation quanti¬ 
ties are defined as gradients of scalars, just as in general 
scalar-vector-tensor decompositions. To simplify nota¬ 
tion and make some calculations easier, without loss of 
generality due to rotational invariance, the authors of [38] 
choose the Fourier wave vectors k to be oriented along 
the z direction. We follow that convention. Therefore, 
these perturbations appear solely in the z direction and 
gradients turn into partial derivatives along z, d^. We 


end up with the perturbations 
tensor: SAj^ = atij6}Sj^ 

vector: = a (Y^, 0,0, dzM^) 

where/, J e {1, 2} 

[ 16 ) 

scalar: SA[ = a(0, SB, 0, 0) 

SAf^ = a(0,0,SB,0) 

SAl = a{dzY, 0,0, SB + djM) 

Sgoo = a^2$G, Sgos = a^d^b 

which show that there are five modes in the scalar sec¬ 
tor, six in the vector sector, and four modes in the tensor 
sector. In the following subsections we will analyze these 
in detail. Not all of those will turn out to be physi¬ 
cally propagating degrees of freedom. Some of them will 
have nondynamical, algebraic, or first-order time deriva¬ 
tive equations of motion and we will refer to those equa¬ 
tions as constraints and nondynamical modes. 

In total, we will see that the theory in Eqn. (1) has 
8 physical degrees of freedom - 2 in the scalar, 2 in the 
vector, and all 4 in the tensor sector. The equations of 
motion for the perturbations follow from the fluid equa¬ 
tion (3). We will treat the three sectors individually in 
the next three subsections. 


A. Tensors 

The tensor perturbations describe gravitational wave 
propagation in this gauge fluid. The behavior of these 
is central to this discussion because the broken parity 
symmetry will become visible in this sector. Left- and 
right handed gravitational waves turn out to obey differ¬ 
ent equations of motion and therefore evolve differently 
from each other. Ultimately this has an imprint on the 
polarization of the CMB. 

To obtain the tensor equations of motion we expand 
the action of Eqn. (I) up to second-order in the per¬ 
turbation quantities and solve the corresponding Euler- 
Lagrange equations. Note that one has to explicitly 
add an external matter Lagrangian to obtain the correct 
equations as gravity couples to all stress energy. Without 
these the Friedmann equations would be sourced solely by 
the gauge field. The tensor perturbations in Eqns. (13), 
tij and hij, are transverse, traceless, synchronous ten¬ 
sors, following Refs. [38, 39]. Additionally, tij and hij 
are invariant under general covariance and SU(2) gauge 
transformations. Therefore, we adapt the usual -I- and x 
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polarization description of gravitational waves with 




/O 0 

0 h+ 

0 hx 

Vo 0 


0 0 \ 
hx 0 I 
-h+ 0 
0 0 / 


which enter the tensor metric perturbations in Eqn. (13). 
Next, we dehne left- and right handed circularly polarized 


h'L + 2-h'„ + 


t'k + 2 —+ 


a^Ml 

tR = - 


k -f — — 
a 


hR 

2 

dMp 


The left-circularly polarized gravitational wave propa¬ 
gation equations are obtained upon the exchange k —>■ 
—k. This has consequences for the evolution — circular 
dichroism — as can be seen by examining the terms with 
a single power of k in the above equations. These terms 
show that the effective mass term —2gYM(f> for the YM 
tensor perturbation t and the coupling between h and t 
differs for left and right circular polarizations. For clar¬ 
ity, we have omitted the anisotropic shear contributed by 
other species, such as photons and neutrinos, although 
these effects are included in our CMB analysis. 

Examining the form of the gravitational wave equa¬ 
tion in the presence of the gauge field, we make the fol¬ 
lowing observations. First, the parity violation depends 
solely on the coupling ^ym and is due to the antisym¬ 
metric Levi-Civita tensor in the internal SU(2) space, 
which we have identified with the principle axes of our 
physical space. Therefore, the coupling constant gyu 
will have a big impact on the chiral properties of this 
model. Second, the tR/p terms in the gravitational wave 
equations may be thought of collectively as represent¬ 
ing an anisotropic shear source. This is similar to the 
sources arising from photons and neutrinos which damp 
subhorizon scale gravitational waves in the standard cos¬ 
mological model. However, here they introduce the par¬ 


wave amplitudes in the standard way. 


hLM = {h+ ± ihx), 


tL,R 


1 

71 


{t+ ±itx), 


upon which the left- and right handed systems separate. 

Following through with the procedure outlined above 
we obtain the tensor perturbation propagation equations: 
the equations of motion for the Fourier amplitudes of 
a right-circularly polarized gravitational wave traveling 
in the +z direction, and the corresponding gauge field 
fluctuation, are given by 


(X 

{k — ) 

[(fc -I- gYM(l>)gYM<P^hR — . (14) 


ity breaking features into cosmology. And finally, we 
note that the background gauge field contributes a novel 
masslike term involving differences in the background 
gauge electric and magnetic fields cx {gYM4’‘^~^'^)- These 
equations will be used later in Sec. IV to describe grav¬ 
itational wave propagation and CMB temperature and 
polarization anisotropy spectra. 

B. Vectors 


2 

aMp 


In this subsection we will derive the vector perturba¬ 
tion equations; we will show that they lead to subdomi¬ 
nant effects in Sec. IV C. Therefore, for the CMB results, 
we will be able to ignore them in our CAMB implemen¬ 
tation. 

The vector perturbations of the metric and the gauge 
field are dehned in Eqn. (13). As it will turn out only 
is a dynamical field. To obtain the equations of motion 
we expand the Lagrangian to second-order in the vec¬ 
tor perturbations, add an external matter part, which 
we choose to be massive scalar, and compute the Euler- 
Lagrange equations for every component of the vector 
perturbation quantities This yields an 

equation for each component of the vectors of which there 
are two. These are 


Yi : k'^iaMiY - ikgYM{(p{aM2y - (j)'aM2) - a{k'^ + 25fYM7)Ti -f 2iakgYM4>Y2 = gYM(t>^{ikB2 - 2gYM4>Bi) 

Bi : MY{6a'^ - a^k‘^)Bi = - 5ym7)^i + 2agYM7(25YM</>Yi - zfcYa) + 2igYMk4>{(t){aM2y - (j)'aM2) 

(15) 

Ml : ik{aMy 2a'M^ a"Mi -|- agYM^Tfi) = iki^oYi)' (7ym(<(’(oT2)^ -I- 2(f>'aY2) — 3ym<('(<(’.S2 T B 2 ) 


and the other set of equations is obtained by assigning {l,2,i} —)• {2,1,—i}. Conveniently the Yj equations can 
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be solved algebraically and eliminated. The Bj equa¬ 
tions, the vector metric perturbations, also couple to any 
other source of vector modes such as the photon fluid. 
For simplicity, we have ignored external sources of vector 
perturbations in deriving the above, middle equation. 


C. Scalars 

The scalar perturbations do not affect the chiral sym¬ 
metry breaking in the theory, but to properly compute 
the impact of this theory on the CMB we need the prop¬ 
agation equations for them. Here we use the definitions 
in Eqn. (13) together with energy conservation of the 
stress energy tensor in Eqn. (4) = 0 and the fluid 

Eqns. (3) to obtain the full set of perturbation equations. 
Defining the transformations 


set of equations that describe the evolution of the fluid 
variables is equivalent to Eqns. (16, 17) and is given by 


a a p 6 3 

6»' = f 7 — - <i)G - cr 


4 p 


(19) 


1 


Sp = ^5p. 


Combining these gives second-order equations again. 
The relationship between pressure and energy density 
perturbations is the same as for radiation. The scalar 
sector of the YM fluid behaves very similarly to regular 
radiation. This can be seen most easily by transform¬ 
ing into the conformal-Newtonian gauge using the usual 
transformation laws 


y = oY, 5m = ak'^M, 5(j) = aSB 

will make these equations simpler. Using energy conser¬ 
vation and the fluid equations, the differential equations 
for the 4 new scalar quantities d4>, y, Sm and b read 

Sm” — k^y' -\- k^dcj) + k'^(j)% — 0 

5(j)'' + k^6(j) - 2gYu^<t)^ {6m - 36(1)) + k^(l)'b (16) 


a' 

o — Wgi — 4^cN 

« ( 20 ) 

b' + —b = 4>gi = xpcN 
a 

where 'kgi and $gi are the gauge-invariant scalar poten¬ 
tials from [44] . The subscript cN stands for the conformal 
Newtonian gauge from [43], which is the gauge that we 
want to translate the perturbations into. With these we 
get 


together with the first-order constraint equations 


5m' - 6(j)' - (fc^ -H 2gYM^<)>^) y + 2gYM^(t>^b - = 0 

(f) {5m — 36(j)) + (py' + 2(j)'y + (/)^4)g — 3(j)'(j)b — (j)‘^b' = 0 

(17) 

The first set of equations (16) describes true physically 
propagating degrees of freedom for 5m and 5(j). The 
second set (Eqns. (17)), however, are merely constraint 
equations that do not describe physical degrees of free¬ 
dom (y, b). They lack second derivatives in a timelike 
component. 

We can also express these in more common fluid vari¬ 
ables to shed some light on these quantities. Equations 
(21) and (22) in [43] define the energy density perturba¬ 
tion 5p, the divergence of the fluid velocity 9, and the 
shear stress a. Reading these off from the perturbed 
stress energy tensor yields 


= -^9 y + 5 ym -^dm - jdm 




a** 

{p+p)0 = 2k‘^^5(j) + 2gYM^k^^y - 2gYM^k^^h 

6p = ^ {36(j)' — 5m') -\- 2yYM^^ {35(5 ~ bm) 

(j)'^ 9 4>' 

( 18 ) 

The first of these equations is the most interesting one 
as it describes the scalar contribution to the shear. The 


{cN) Y = -^0 + 4</>'^ 

{cN) 9' = ^k'^6 + k^ip — k^a, 


( 21 ) 


where the notation “cN” is to remind us that all variables 
are now in the conformal-Newtonian gauge. This is iden¬ 
tical to Eqns. (64) in [43]. It is just the contribution to a 
in Eqn. (18) that gives deviations from normal radiation 
in the scalar sector. 

In the case yYM = 0, the evolution of the scalar shear, 

{cN) a' = ^9, (22) 

may be obtained from the perturbation equations. In 
the case yvM 0, special care must be given to evolve 
the scalar shear as will be explained in Appendix A. To 
implement the code in GAME we follow the definitions 
of the perfect fluid perturbations from [43] in the syn¬ 
chronous gauge. We keep using adiabatic initial condi¬ 
tions for the density contrast <5 upon introducing the YM 
fluid. 


IV. SOLUTIONS OF THE PERTURBATION 
EQUATIONS 

With the evolution equations derived in the previous 
section, we can analyze the modified behavior of grav¬ 
itational waves as well as their impact on the CMB. In 


7 


the following sections we will study the cosmological evo¬ 
lution of gravitational waves and present a technique to 
efficiently obtain those results. Finally, we will determine 
the impact the YM fluid has on the power spectra. 


A. Prologue: Secular perturbations 

In the following section we will derive the scalar per¬ 
turbation equations and compare these to the standard 
ACDM results in the conformal-Newtonian gauge. Be¬ 
fore we make the formal derivation we will briefly review 
secular instabilities as an example how perturbations can 
grow indefinitely without actually affecting the physics. 

Consider a pendulum of length I swinging at an angle 
0 relative to the vertical. The full ordinary differential 
equation describing this system is 

0 + At sin 0 = 0 (23) 

where At = g/l. Note that expanding this equation up to 
third order in 0 yields an equation similar to the equation 
of motion for our gauge scalar (see Eqn. (8)). Introducing 
a perturbation 0 = 0o + 60 yields 

60 + K cos 0o60 = 0 , (24) 

where 0o follows the background trajectory. This is the 
linearized perturbation equation for this system. Now, as 
cos 00 oscillates, the perturbation 60 grows linearly which 
looks like an instability of the theory. This is tradition¬ 
ally called a secular perturbation. However, obviously 
the simple pendulum does not contain any unstable fluc¬ 
tuation. The actual problem lies in misinterpreting initial 
conditions for this perturbation [40]. Shifting the initial 
conditions of the pendulum propagation might look like 
this growing mode. Secular perturbations often refer to 
these shifts in initial conditions and are therefore not a 
true physical instability. 

To be explicit, let us consider a small perturbation to 
the gauge field solution (9) , 0 -I- 6(j), and linearize 

the equation of motion (8), whereby 

6(f>'' + 6gYyi(j3^6(j) = 0. (25) 

The solution to the above equation is 

6(j) = (^6(j){Ti)(l)' + ^6(j)'{Ti){T(j)' + (())^ . (26) 

Although (5(/) contains a term that grows linearly in t, it is 
obvious that this signals no physical instability; instead, 
it indicates a limit to our linearization approximation. If 
we simply adjust the initial conditions for (j) or else refrain 
from linearizing the perturbation equation, then we can 
make 6(j) well behaved. 

We will encounter these explicitly in the solutions to 
the vector perturbations in Sec. IV C and scalar pertur¬ 
bations in Sec. V. 


B. Gravitational Waves 

In this chapter we will use the formal results from 
Sec. Ill A to describe how gravitational waves behave dif¬ 
ferently with the YM fluid present compared to standard 
ACDM. Most prominently the fluid introduces chiral ef¬ 
fects into the tensor sector of gravity and therefore affects 
the behavior of left- vs right-handed polarized gravita¬ 
tional waves. This cosmological parity rotation is man¬ 
ifest in the coupling to gravitational waves. To visual¬ 
ize this, consider a gravitational wave passing through 
the gauge field described in the above sections. Gen¬ 
erally speaking, gravitational waves induce a quadrupo- 
lar distortion, alternately squeezing and stretching the 
stress and energy of the YM field. However, the gauge 
field itself possesses a preferred handedness via the right 
handed SU(2) structure constants. Since the fields are 
flavor-space locked, where the principal axes of SU(2) are 
aligned with the spatial coordinate basis, the gauge field 
stress energy will vibrate in sympathy to a right handed 
wave, and with antipathy to a left handed wave. A nice 
visualization for a similar effect is the rattleback top [45] . 

We now analyze the evolution of gravitational waves 
described in Sec. HI. There is a rich variety of behav¬ 
ior in the evolution of this system, dependent upon the 
coupling gvM, the abundance Rym = PYivi/prad during 
the radiation era, the relative contributions of electric 
and magnetic field energy, (jJ and and the initial 

conditions for the perturbations ^r/l- For these 

purposes, we omit the anisotropic shear contributed by 
other species, such as photons and neutrinos, although 
these effects are included in our CMB analysis. Further¬ 
more we assume a standard scale-free primordial tensor 
power spectrum. 


1. Long wavelengths 

We first tackle the long wavelength behavior of the 
gravitational waves. For simplicity, we begin our investi¬ 
gation with the simpler scenario, setting ^ym = 0 which 
corresponds to electrodynamics with three flavors or col¬ 
ors, also known as “color electrodynamics”. In this case, 
the distinction between L— and i?— handed gravitational 
waves is gone, so we label the tensor perturbations with 
the subscript “A” for ambidextrous. The equations for 
Ha and Ia (14) may be combined as 

h'X + 2°^h'A + ik‘^-2kl^)hA 
= —4A:q^ [ dr'h A cos k{T — t'). (27) 

“ dri 

Hence, there is a time-dependent effective mass term, 
miff = — 2fcoao/a^. In the case k < y^koaoja this 

term is negative, producing an enhancement of the grav¬ 
itational wave amplitude. To see this, we take the long 
wavelength limit, fc —> 0, in which case the equation be- 



comes 


(28) 


and inserting it into the first equation, then we find at 
the fixed point 


/i'; + 2-h(4 + 2kmhA = 


where hA{Ti) is the initial value of the wave ampli¬ 
tude. By inspection we note that there is a fixed point 
where the amplitude approaches a constant value, Ha —t 
2hA(Ti). We show this more rigorously by assuming a 
radiation-dominated expansion scale factor a = airjTi 
and, assuming the standard initial conditions = 0, 

then 


hA{T) 


hAin) 


2 -|- Cl 



ni,2 = - ^ (^1 ± \/l - SiaokoTi/ai)"^^ 


(29) 

(30) 


where Ci = —(l+ni)/(l + 2ni). Since the exponents ni ,2 
are both negative, the time-dependent terms decay, and 
the fixed point is soon reached. It is clear that for modes 
with wave number fcro ^ 1 and also k '/^koao/a, 
the amplitude of long wavelength gravitational waves is 
doubled. 

We can extend the analysis of the behavior of long 
wavelength gravitational waves to the coupled, YM case. 
To begin, we set /c —)■ 0, and assume a pure radiation 
expansion rate so that the a"/a terms may be neglected. 
Next, we introduce a change of variable u = at and in¬ 
troduce the new time variable x = kg{T — Ti) + C 2 - The 
system of equations (14) now becomes 


cPh 

dx"^ 


d?u 

dx'^ 


^dlnadh 2 . , ,on . 

-^ (V' - ) h 


dx dx 
2 

a^A 


i2A2 


2 

W 


, , du 1 dl^ip \ 

+ 2 

(31) 

dh 1 d'^ip \ 
dx 2 dx"^ ) 

(32) 


where A = gYMO,oMp/kg. The functions ip and dip are the 
background field and derivative scaled to unit amplitude 
as 


iP 

dip 


5ym 


kg 

gvM 

-kf 


(p = sn(a;| — 1) 

(p' = cn(a:| — 1) dn(x| — 1). 


(33) 


In the case that ip, dip are slowly evolving, effectively con¬ 
stant, then we note there is a self-consistent fixed point 
solution whereby h ^ hf and u ^ atf such that h'j and 
are negligible. From Eqn. (31) we determine 


dip'^^=-\{iP^-diP^)hf. (34) 

Next, integrating Eqn. (32) so that 



(36) 

These results are confirmed by a numerical integration. 
To reformulate these results in terms of the initial condi¬ 
tions for Ip and dip, we obtain 

h{Tf) = 2 siiP^ 9 h{Ti) 

t{Tf) = —2(1 -I- Zf)koTf sin 9 cos 29 h{Ti) 

where Tf,Zf are the conformal time and redshift and 9 is 
defined in Eqn. (11). These assume standard initial con¬ 
ditions for h and we have set the initial perturbation in 
the YM fluid t{Ti) = 0. In the case 9 = 7r/4, the effective 
mass term in the evolution equation for h, proportional 
to the difference between the background field electric 
and magnetic energy densities, vanishes, so that h and t 
do not evolve on long wavelengths. This is the standard 
case, which we focus on in this paper. 

In the case 9 = ti 12, the ip^ term drops out of the effec¬ 
tive mass term, just as in the case of color electrodynam¬ 
ics where gYM = 0, and the amplitude of h doubles. The 
most surprising case, however, seems to be when 9 — 0, 
so that the effective mass is entirely due to the ip'^ term, 
which acts to damp the gravitational wave amplitude. 
Assuming that the evolution begins at some time shortly 
after inflation, then the fixed point is soon reached well 
in advance of horizon entry by any wavelengths of inter¬ 
est. The prediction, borne out by numerical integration, 
is that the gravitational wave and tensor modes of the 
gauge field are damped out. (We note that such suppres¬ 
sion could modify the predictions of the gauge-flation and 
chromo-natural inflation scenarios [36, 37, 39] or other 
inflationary models that traditionally overproduce grav¬ 
itational waves. This effect has not been explored before 
and might mitigate the shortcomings of these theories 
that ultimately led to the predictions that ruled them 
out.) 

This solution (Eqn. (37)) also allows us to predict the 
amplitude of superhorizon modes at some time late in the 
radiation-dominated era, which we may then use as the 
initial data for a numerical study as the modes proceed 
to enter the horizon. 

We now focus on a minimal scenario in which the ini¬ 
tial field energy of the YM fluid is split equally between 
the electric and magnetic field, (p' = gYM<P'^, whereby 
9 = tt/A. These initial conditions live closest to the stan¬ 
dard cosmological model, as long wavelength modes are 
unaffected by the gauge field. We further assume equal 
amplitude scale-free primordial spectra of left- and right 
handed gravitational waves. These initial conditions, and 
the assumption that the initial tensor fluctuations of the 
gauge field vanish deep in the radiation era, allows the 
YM fluid to behave like radiation at early times. The 
effects of the gauge field on the subsequent evolution of 
the gravitational waves are illustrated in Figs 1 and 2. 
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10-2 102 10 "^ 


conformal time: r 

FIG. 1. Gravitational wave amplitude evolution as a function 
of conformal time is shown for the case ^ym = 0, Rym = 
0.03 (blue) and wave number fe = 10 h/Mpc, as compared 
to the standard case 7?ym = 0 (green). The excitations of 
the gauge field are shown (red) as an offset minus a constant 
times (atA)^, to illustrate their complementary behavior. The 
solid (black) lines show the results of WKB solutions for the 
envelopes of the oscillatory waveforms. 



10-2 10 *^ 102 104 

conformal time: t 


FIG. 2. As in Fig. 1, the gravitational wave amplitude evo¬ 
lution as a function of conformal time is shown, this time for 
the case qym = 10“®°, Rym = 0.03 (blue) and wave number 
k = IQ h/Mpc. The waveforms capped with solid lines are 
right handed; those with dashed lines are left handed. 


The evolution of the gravitational wave amplitude is 
shown for a variety of cases in Figs. 1 and 2. We begin by 
examining the behavior in the case qym = 0 , correspond¬ 
ing to color electrodynamics. The background solution 
has (j)' constant, so the Ha evolution equation (where the 
subscript “A” is for ambidextrous, since there is no par¬ 
ity violation in this case) has a tachyonic mass that is 
responsible for the growth of long wavelength modes. As 
we have shown analytically, Ha doubles for modes out¬ 


side the horizon, relative to the standard case. As seen 
in Fig. 1, the amplitude of (a/i^)^ (blue) is 2^ times the 
standard case (a/i)^ (green) going in to the first oscilla¬ 
tion. We also notice that after modes enter the horizon, 
there is a slow exchange of amplitude between h and t. 


2. Short wavelengths 


To investigate further, we make a Wentzel-Kramers- 
Brillouin (WKB) analysis [46] to obtain the gravitational 
wave behavior up to high wave numbers - a configura¬ 
tion space region that is unaccessible numerically due to 
high computation cost. The behavior of the gravitational 
wave amplitude as a function of time for different values 
of the wave number may be obtained by solving Eqn. (27) 
numerically. However, once a mode enters the horizon, 
with kr ^ 1, then the speed of computation slows. In 
this case, we may adopt the WKB approximation devel¬ 
oped here in order to solve for the more slowly evolving 
envelope of the wave amplitude. We write 


Ha = ^ sin(fcr -I- 0^), Ia = ^ ^ sin(fcT -I- 0^) (38) 

and define A9 = drj — 9u- The functions we have to 
determine are rj, v and AO. In the limit fcr ^ 1 the 
equations of motion become 


r]'=—kfj — vcos AO, v'= k^—rj cos AO 
a a 

AO' = ko— (---] sinA 6 ». 


a 


V 


(39) 


Hence, our procedure is to solve the full set of equations 
until fcr ^ 1 , at which point we set initial conditions for 
the envelopes 77 , v and the phase separation AO. This 
set of equations is more easily solved and provides an 
excellent fit to the full numerical solution. First of all, 
we notice that the sum 77 - 1 - 7 ^ is a constant, indicative of a 
conserved quantity in the coupled system of gravitational 
waves and gauge field tensor fluctuations. This visualizes 
the property of Fig. 1 nicely: changes in the envelope 
of hA are compensated by the envelope of Ia- We plot 
the two solutions together in order to clearly show this 
effect (see Fig. 1). Second, we notice that each time 
derivative term contains a factor koag/a, suggesting that 
the natural time parameter for the modulation of the 
envelope is k^dra^f/a. Our solutions show agreement with 
the numerical results, with the amplitude of Ha and Ia 
changing with a period defined by fco f draoja = 27 rn for 
77 , an integer. 

The high-frequency behavior of the gravitational waves 
in the case ^ym 7 ^ 0 is also amenable to a WKB analysis. 
Upon making the same definitions for 77 , ly, the equations 
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of motion become 

r]'ji = -kQ^VR^dip cos Mr- iip- sin Mr ) 

v'r = ko°^r]R{dipcosA6R - sinA0ij) 

+ ko— (— - 

a \Vb. JJrJ 

X (di/'sin cos A0i{) 

ti't = — fcg —J^ifdV’cos A0 l + sin A0i) 
a 

v'r = ko —7]L{dilj cos AOr + ij}'^ sin A0 l) 

A0^ = -fcgV' + fco- 

a \vl vrJ 

X {d-ip sin A9r — ip'^ cos A9r) 

Our procedure is to solve the full set of equations, for 
both R— and L—handed gravitational waves, until fcr 3> 
1. Then we use our numerical solution to set initial condi¬ 
tions for 77, V , A9 and then evolve the envelopes forward 
to the present day. The results are qualitatively similar to 
the ambidextrous case, except that for certain ranges in 
the value of the wave number k there is a large difference 
in the amplitudes for R— and L—handed gravitational 
waves. Looking at the figure showing the evolution of 
the wave amplitude as a function of time (Fig. 2), we 
again see a slow mixing of power between the gravita¬ 
tional wave and the tensor fluctuation of the gauge field. 
There is a slight boost to the i?—handed gravitational 
wave at horizon entry relative to the standard case, and 
a small suppression of the L—handed wave. The WKB 
approximation is shown to do an excellent job of tracing 
the shape of the envelope in the fast oscillating regime. 

3. Gravitational wave spectral density 

The gravitational wave spectral density is the gravita¬ 
tional wave energy density per log frequency interval, in 
units of the critical density. Following notation in [12], 
the spectral density for inflationary gravitational waves 
is 

s—L,R 

where the inflationary initial conditions for the amplitude 
are Aj^ ^ i / Mp)^ / tt and T is the transfer function 

defined such that T(r, k) = h{T, k)/h{Ti, k). The gravita¬ 
tional wave spectral density flow is shown in Figs. 3 and 
4, where the amplification and periodic modulation are 
clearly seen. The WKB solution predicts a peak or dip in 
the spectral density every 5 orders of magnitude in k for 
Rym = 0.03; the origin of this oscillation in the spectral 
density is the same as the oscillation of the amplitude of 
the envelope of the wave amplitude between the metric 
and gauge field tensor perturbations. 


The difference in the evolution for left- and right- 
circularly polarized waves is primarily due to the effec¬ 
tive mass term for the gauge field tensor perturbations, 
—2kgYM(p, which is tachyonic for right handed modes. 
This behavior leads to an interesting effect: The growth 
(suppression) of tp {tp) is transferred to hp (hp) as the 
mode enters the horizon (see Fig. 2). Once the rel¬ 
ative amplitude is locked in at horizon entry and the 
fields begin to oscillate rapidly, the slow exchange of am¬ 
plitude between hp/p and tp/p again comes into play. 
The WKB analysis for subhorizon modes shows that 
~k^‘R/L oc 1/a^ and the exchange is oscillatory with 
similar phase if the background field is not yet oscillatory. 



10"^ 10“ 10' 10" 10'’ 10** 

k [inv-Mpc] 


FIG. 3. The gravitational wave energy density spectrnm 
is shown as a function of the comoving wave number. An 
ambidextrous, scale free spectrum at an inflationary scale 
Hi — lQ~^Mp is assumed. The present-day spectrum in 
the case with ^ym = 0, J?ym = 0.03 displays large oscilla¬ 
tory features due to the coupling between the gravitational 
waves and the gauge field. For comparison, the standard case 
without the gauge field is shown, as well as the effect of Stan¬ 
dard Model particle free streaming and freeze-outs (repro¬ 
duced from Ref. [12]). 

The WKB solution allows us to easily calculate the 
gravitational wave spectral density out to high frequen¬ 
cies. Looking at Fig. 3, and comparing with the standard 
case obtained from [12] we notice several significant fea¬ 
tures. First, the long wavelength amplitude is higher 
by a factor of 4 than in the standard case, as expected. 
Second, there are slow, secular oscillations in the ampli¬ 
tude, as predicted based on the envelope oscillation phase 
aofco / dr ja. The decay of the oscillation amplitude at 
high-frequency is due to the choice of starting time for 
our integrations. Those wave modes start outside the 
horizon, but did not reach the amplitude-doubling fixed 
point before horizon entry. We started these integrations 
at a redshift z ~ 10^® for reasons of numerical accuracy. 
If we had started at a more physically realistic z ~ 10^^, 
then the decay of oscillations in ^aw would be shifted to 
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FIG. 4. As in Fig. 3, the gravitational wave energy den¬ 
sity spectrum is shown as a function of the comoving wave 
number, but with (Jym = 10“®°, Rym = 0.03. The long 
wavelength modes are unaffected given our choice of initial 
conditions in the background field. However, there is a strik¬ 
ing difference in the spectra of L- and i?-handed gravitational 
waves. 


even higher frequencies. Third, and finally, for the values 
of the parameters ^ym and Rym used in this figure, the 
oscillations in the amplitude of flow due to the coupling 
with the gauge field are much larger than the imprint of 
the thermal history of the relativistic fluid as calculated 
in [12]. 

The gravitational wave spectral density for the case 
of the fully nonlinear Yang-Mills gauge field is shown 
in Fig. 4. Here we see the striking asymmetry in the 
left- and right handed gravitational wave spectrum near 
k ~ 10°, 10° inv-Mpc. Although our numerical calcula¬ 
tion extends only to fc = 10° inv-Mpc, the asymmetric 
modulation of the two spectra should continue to higher 
wave numbers. Should the direct detection of a stochastic 
gravitational wave background ever threaten to become 
a reality, one might consider attempting to distinguish 
between the left- and right handed spectrum. Finally, 
we reiterate that the details of our spectrum depend on 
choices of initial conditions for the tensor modes 
(which we have set to zero) and the background field (j) 
(in terms of the split between electric and magnetic en¬ 
ergy). 


C. Vector modes 

Here we examine the behavior of the vector modes, 
based on Eqns. (15). We solve for the constraint modes 
Yj as described earlier, leaving second-order equations 
for Mj that are sourced by Bj. 

For our first attempt to study this complicated sys¬ 
tem of equations, we set Bj = 0. The Mj differential 


equations are coupled, so we make the transformation 
Ml = 2gYyik‘^a‘^fi + igyyika'^ {k^ + ^ ( 42 ) 


to unmix the second-order terms, where // are the new 
vector perturbation variables and again, the mapping 
{1, 2, i} —>• {2,1, —i} yields the second equation. The new 
equations are still rather complicated, so we consider the 
limiting cases of very high and very low frequencies. In 
the case k ^ 5ym 0, the equation of motion becomes 

fi + (i-fi + k^fi = 0 (43) 

a 

which describes a damped harmonic oscillator with the 
damping depending on the cosmology. In the case of 
power-law expansion, with a oc where w is 

the background equation of state, then // oc jnikT)/T^ 
and n = (5 — 3w)/(l -|- 3w). To estimate the conse¬ 
quences of this solution, one would typically compute 
the density contrast Sp/p or a, similar quantity to see if it 
grows or decays. However, in this specific gauge choice of 
Eqn. (13) the energy density and pressure perturbations 
only appear at second-order. The off-diagonal terms in 
the YM stress energy tensor are first-order in perturba¬ 
tions though, so we compare the off-diagonal terms to the 
unperturbed energy density to determine if the pertur¬ 
bations grow. Eor example, the dominant contributor to 
the momentum density is 6Ttx —igYMkM 2 (j)(l)'/a and 
Ml ~ igYM^‘^k^h/4’- Using our analytic solution, and 
Ttt = a^Prj} then 


STtx 

Ttt 


2gYM^^^' 
3((/)'2 -b 


a^T '^jr 


(kr) 


Akz 


(44) 


Using the asymptotic behavior of the spherical Bessel 
function, jn{x) —>■ 4 cos (a; — {n + l)7r/2) as a; —)■ oo, 
then we see that the momentum density conveyed by 
the vector perturbation oscillates with constant profile, 
^Ttx/Ttt oc Similar results are obtained for other 

nonzero off-diagonal components of the stress energy ten¬ 
sor. Hence, there is no gravitational instability for modes 
with k ^ gYWp- Based on typical numbers used in our 
study, this translates to (/ym')' I0“°iFo- For practical 
purposes, this applies to all modes of interest. 

Next we include the vector metric perturbations, by 
solving the algebraic constraints for Bj and Y/, leaving 
second-order equations for Mj. The same transformation 
is made to unmix the the second-order terms in favor of 
the variables fj. In this case, the high-frequency limit for 
the evolution of // is the same as before. Hence, our anal¬ 
ysis from above still holds true, that the high-frequency 
vector perturbations oscillate with constant profile. 

To study the very low frequency vector perturbations, 
we start again with the case Bj — 0 and proceed to unmix 
the equations of motion. Taking the leading terms in 
the limit k gYM</>) then the equation of motion for // 
reduces to {Atpfi)" = 0. Since the dominant contributor 
to the momentum density perturbation in this limit is 














12 


STtx ^ igYyik^{a^<j)‘^fiy/a'^, then in units of the energy 
density we find 



3(0'2 + 


+ C2(r(^' + '/')), 


(45) 


where ci ^2 are integration constants. This growing term 
appears to indicate an instability due to the linear growth 
in r; however, it is simply the secular perturbation dis¬ 
cussed in the prologue to this section. That is, by ad¬ 
justing the initial values of (j) and (j)' then the constants 
Cl .2 can be made to vanish. When we include the vec¬ 
tor metric perturbations, once again the same equations 
of motion and momentum density perturbation are ob¬ 
tained in the small-/c limit. There is no instability. 

Because there is no external source of vector pertur¬ 
bations, and because the YM vector does not amplify 
any perturbations, we have chosen to omit computing 
the contributions of vector perturbations to the CMB 
anisotropy in Sec. V. 


D. Scalar perturbations 

The scalar perturbations are stable, as we now illus¬ 
trate. The equations of motion for the scalar degrees of 
freedom, 5(j) and dm, can be simplified by making the 
change of variables 


= Ai, dm = Ai - A2y^P~+2g^J^/(gYM0)- (46) 

In this case, the second-order terms in the equations of 
motion unmix, so that 

A" + MyAj=S, (47) 


where 


Aiii — -|- 45 ^ 1 ^ 1 *?^^ 

Mi 2 = M 21 = 2gYM4'\J k? + 2g‘^y^(jP 

M22 = k^ + 2g^YM^^+ 


Si = 4g^M<(>^4’G - ^'^>0 “ 

^ (3A;^ -b 2gl^(j)‘^)gYM(l>(t>' _ 9Yyi(t)^ 


(fc2 + 2g^j^</)2)3/2 

(fc^ -b Ag‘i^yi4P‘)gYM(t>^ 
^/kP- + 2gYM4>^ 

ffYM'/'V 


k'^ + 25ym'^^ 
2k‘^ gYM(t>''^ ^ 


k'^b' 


(P + 2g|^<(.2)3/2 


4>r 






These are just two coupled, driven harmonic oscillators. 
The eigenvalues of the matrix M are positive definite, so 
the oscillator is stable. At long wavelengths, however, we 
can see the imprint of the secular instability. In the fc —>■ 0 
limit, the homogeneous equations become 6m" = 0 and 


6u" + 6gYM'P'^^'"' = 0 where Su = 3S4> — Sm. Of course, we 
recognize the Su equation as Eqn. (25). When we insert 
the growing solution (26) into, say, the energy density in 
Eqn. (18) then we find Sp/p oc 6u'{Ti) is a constant and 
ought to be absorbed into the background. Similarly, the 
shear appears to grow linearly 

^9YM(l^^i^rn{n) + Sm'{Ti)T) - (ji'Sm'in) 

but again this is an artifact of our linearization. The 
system is stable as shown by our numerical calculations 
in the next section. 


V. CMB 

Ultimately, detecting the imprint of this chiral YM 
fluid requires an analysis of the CMB sky. In this section 
we compute the power spectra and describe ways to de¬ 
tect deviations from vanilla ACDM. To evaluate the im¬ 
pact of this scenario on the CMB, we have implemented 
the scalar and tensor perturbations of the gauge field 
into CAMB [30]. We implement the standard adiabatic 
perturbation for scalar and tensor perturbations but set 
t{Ti) — 0 lacking an early Universe theory for this toy 
model. This leaves us with solutions that live close to 
vanilla ACDM so these new effects are not artificially en¬ 
hanced. The gauge field has the biggest impact on ten¬ 
sor correlations. The scalar sector also receives correc¬ 
tions due to the gauge field, in the form of an anisotropic 
scalar shear, but the impact on the scalar CMB spec¬ 
trum is small. We ignore the vector perturbations which 
we showed to decay in Sec. IIIB. We use the parameter 
Rym describing the ratio between the YM fluid density 
and the total relativistic energy density. We assume the 
fraction of critical density in the relativistic fluid is fixed 
by slightly adjusting the sum of the neutrino masses upon 
introducing the gauge field. This avoids simply compar¬ 
ing the YM fluid in terms of an increased effective neu¬ 
trino number N^s as we are trading the sum of total 
neutrino mass for ANgg. We otherwise assume standard 
ACDM parameters. Details of our implementation can 
be found in Appendix A. 

The CMB polarization can be decomposed into gra¬ 
dient U-modes and curl U-modes. In the tensor sector 
the gauge field introduces two main effects. First, the 
left- and right handed contributions to the BB spectrum 
now differ, as shown in Fig. 5. Hence, the temperature 
and polarization anisotropy due to gravitational waves on 
roughly degree scales is dominated by a superposition of 
right-circularly polarized gravitational waves, which im¬ 
print left-helical patterns on the sky. It is curious to see 
that the individual contributions deviate strongly from 
vanilla ACDM but conspire in a way that puts the com¬ 
bination of both close to the expected standard cosmol¬ 
ogy result. This holds true for a wide parameter range, 
and these effects can be seen in Fig. 6, the correspond¬ 
ing graph with ^ym = 10“®°. The general features in 
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FIG. 5. CMB BB polarization autocorrelation spectra. The 
pure left- and right handed contributions deviate strongly 
from ACDM (dashed) while their sum is closer to it. Solid 
lines include a gauge field with qym = 10“®®, Rym = 0.1 and 
tensor-to-scalar ratio of r = 0.1. The thin line shows the BB 
spectrum for the case qym = 4 x 10“®®. 
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FIG. 6. Graph analogous to Fig. 5: the gauge field now has 
Qym ~ 10“®®, 7?ym =0.1 and tensor-to-scalar ratio of r = 
0.1. Again, the thin line shows the BB spectrum for the case 
Qym = 4 X 10“®®. 


a smoking gun for chiral effects in the Universe. Typ¬ 
ical predictions of our model are plotted in Fig. 7. In 
that plot we include all types of correlations to make a 
comparison between the signal strength of the different 
spectra easier. ACDM spectra are plotted with dashed 
lines. 

So far this discussion was centered around choosing 
initial conditions that resemble vanilla ACDM. In par¬ 
ticular, the angle 9 in Eqns. (11, 37) was chosen to be 
equal to 7r/4. Deviating from this has two effects: The 
background evolution of (p changes drastically and the 
effective mass term in the evolution equation for h no 
longer vanishes as described in Sec. IVBl. This behav- 



FIG. 7. CMB temperature and polarization auto- and cross¬ 
correlation power spectra. The parity violation allows for T B 
and EB cross-correlations. The dashed lines represent the 
standard ACDM cosmology, solid lines include a gauge field 
with Qym = 10“®®, Rym = 0.1, 9 — 7r/4 and tensor-to-scalar 
ratio of r = 0.1. The TB and EB cross-correlations only 
appear when the gauge field is present. 


this parametrization are similar, although the YM modi¬ 
fied spectrum follows the vanilla ACDM one more closely, 
especially at large scales. Second, because temperature 
T and gradient polarization E are both parity even but 
curl polarization is parity odd, the correlation spectra 
between these types vanish in the vanilla ACDM sce¬ 
nario. However, the parity violation introduced by the 
YM fluid allows for correlations between these types, TB 
and EB [15]. Detecting these exotic cross-correlations is 


ior is shown in Fig. 8. The solid black line describes 
the parametrization that has been employed in the rest 
of this paper, using 9 = 7r/4. Going to larger angles 
9 ^ 1112 lifts the entire BB power spectrum almost by 
an entire order of magnitude. Half of the logarithmic 
increase comes from the new initial conditions, as evi¬ 
denced by Eqns. (37). We are using this equation to set 
up initial conditions for CAMB at a = 10“®. The other 
half is caused by the change in the background field be¬ 
havior. Propagation of the tensor perturbations through 
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this altered background field yields this strong effect on 
the spectra, decreasing the angle lowers the spectrum 
indefinitely. This analysis shows the extent to which 0 
changes the behavior of the YM fluid and how sticking 
to 0 — tt/A gives results closest to vanilla ACDM. How¬ 
ever, to explore the entire parameter range, we need to 
take variations in the distribution of initial electric and 
magnetic energy contribution into consideration. 
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FIG. 8. Curl polarization power spectra. The cosmological 
parameters as well as {gvM, Aym} = {10“®®,0.1} are kept 
constant. The angle 6 is being varied which changes the con- 
tribntion of initial electric and magnetic YM energy. The 
6 — tv/ A curve is closest to ACDM. 

There are many challenges to detecting the parity- 
violating cross-correlations, not to mention the primor¬ 
dial d-mode signal. Galactic foregrounds, magnetic 
fields, weak leasing, and other systematic effects can all 
produce a false positive; fortunately there is no funda¬ 
mental barrier that would prevent a detection that can 
distinguish a primordial signal. (See Ref. [47] for a recent 
summary.) 

But there are other, competing, phenomena that could 
produce a parity-violating signal. First of all, cosmolog¬ 
ical birefringence (CB) can lead to TB and EB power 
spectra by rotating E into B through a novel coupling be¬ 
tween electromagnetism and a cosmic pseudoscalar such 
as quintessence [14]. A second possibility, broadly char¬ 
acterized as chiral gravity, posits a modification of gravity 
whereby an asymmetry between left- and right-circularly 
polarized waves is imprinted on the primordial spectrum. 
The third possibility, as we have shown, is essentially cos¬ 
mic circular dichroism, whereby the asymmetry develops 
with time from an initially symmetric primordial spec¬ 
trum. 

Would an actual detection directly point to chiral sym¬ 
metry breaking on cosmological scales? Do the birefrin¬ 
gence effects mask true chiral physics? In Ref. [48] it 
was shown that the TB and EB spectra can be used 
to distinguish these CB effects from chiral physics. As 
CB rotates the E into a B contribution the measured 


B spectrum would resemble the E one which makes this 
separation into CB and chirality effects feasible. In turn, 
putting limits on the amplitude of these spectra will put 
constraints on chiral physics in general and our model in 
particular. 

In what follows we compute the constraints that cur¬ 
rent and future CMB experiments would put on the pa¬ 
rameters of our model under the assumption that TB 
and EB cross-correlations are measured. 


VI. FISHER FORECASTS 


The deviations seen in the CMB spectra would clearly 
have an impact on the interpretation of a precision mea¬ 
surement of B modes [1]. As can be seen from Fig. 5, 
the gauge field can vary the height of the BB spectrum 
at the reionization bump near i ^ 10 and at the pri¬ 
mary acoustic peak near £ ~ 100 by as much as ±50%. 
However, we have the greatest leverage on new physics 
by focusing on the exotic cross-correlations. Hence, we 
forecast the parameter constraints (T_Rym '^sym using 
these spectra. 

The Fisher-matrix technique gives an estimate for the 
statistical error in a given measurement. Therefore, sys¬ 
tematic deviations are solely contained in the choice of 
a fiducial model along with any theoretical bias effects. 
The Fisher matrix for this case reads 


T 

2^2^ dp dp H JxF 

I X,Y 


(49) 


where p = (RyMjS'ym) ± Pmsmo and X,Y = 
{TB,EB} and S is the Ci covariance matrix. Here, 
the eight cosmological parameters are Pcosmo = 
(uJb,u!c,uj,j,flK,Ho,w,nt,r). (In CAMB, this corre¬ 
sponds to tensor_parameterization = 1.) The Fisher 
matrix E is the inverse of the covariance matrix be¬ 
tween i?YM and ^ym- The derivatives of the are 
obtained using CAMB. We center the derivatives around 
the following fiducial model: We choose ^ym = 10“®®, 
Rym = 0.1, and the standard Planck ACDM values for 
the cosmological parameters [49]. Furthermore, we set 
the angle that distributes initial energy contributions 
between electric and magnetic parts of the YM fluid 
9 = tt/A (Eqn. (11)). 

The matrix is the inverse of the TB-EB co- 

variance matrix given by 




^XiX^ ^X2X4 _|_ ^XiX4 ^X2X3 
21 + 1 


where 


and X = {T, E, B} [48], where W is the window func¬ 
tion and w describes the instrumental noise. Note the 
usage of the superscripts: Here X refers to one type of 
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perturbation only, whereas above it describes the cross¬ 
correlations. The properties of the detector are imper¬ 
ative in the determination of the parameter constraints. 
The instrumental parameters enter the window function 
in two ways: due to the beam width via 

Wl> exp {-l^al/2) 
and the instrumental noise where 

= Aira^/Npix 

and Wgg = = dTrcTp/iVpix with the cross-correlation 

contributions vanishing as the noise in the polarization 
is assumed to have no correlation to the noise in the 
temperature. 

In the window function ai , = 0fwhm/'s/8 In 2 where 
the beam width is measured in radians. Similarly, the 
number of pixels is iVpix = and ap 

are the temperature and polarization pixel noise. These 
are given by = (NET)^A^pix/tobs and ap = \f2aT 
with NET being the noise-equivalent temperature and 
tobs being the observation time. 

The parameters for this analysis are taken from 
Ref. [48, 50] and summarized in Table I. 


Instrument 

^FWHM [arcmin] NET [^Ky^ 

tobs [y] 

Planck 

7.1 

45 

2 

CV limited 

5 

0 

1.2 


TABLE I. Instrumental parameters for the two experiments 
considered in this paper. The parameters are the beamwidth 
^FWHM, noise-equivalent temperature NET, and observation 
time fobs. 


The ID marginalized confidence limits in a scenario in 
which the Planck satellite measures TB and EB correla¬ 
tions are (TgYM = 9.5 x 10“®^ and ap,^^ = 0.030. Eor the 
cosmic variance (CV) limited experiment these numbers 
reduce to crgY„ = 3.4 x 10“®^ and ap^„ = 8.1 x 10“^, 
which would enable us to make a nonzero detection of 
5ym- The 1- and 2-a contours are plotted in Eig. 9 where 
we marginalize over the other parameters. Eor this fidu¬ 
cial model Planck could make a 2-cr detection of Rym, 
but cannot exclude the ambidextrous case since gym = 0 
lies within its 1-a contour. The future looks brighter for 
proposed satellite missions that gets closer to a CV lim¬ 
ited experiment (Prism [51], CMBPol [52]), which could 
put constraints on the chiral asymmetry; for the hducial 
model, the coupling ^ym could be distinguished from zero 
at better than the 2-cr level. 


VII. CONCLUSIONS AND OUTLOOK 

We present a simple model of a dark fluid that breaks 
chiral parity on cosmological scales. We illustrate the 
impact of this model on the gravitational wave spec¬ 
trum and the CMB, computing the power spectra along 



FIG. 9. Forecasted 1- and 2-g C.L. contours under the con¬ 
dition that TB and EB cross-correlations are detected by 
the respective experiments (indicated by the asterisk *). The 
fiducial model is indicated by the black dot which represents 
flYM = 10 and Rym = 0.1 


with new TB and EB correlations that emerge in parity¬ 
breaking models. For a wide range of parameters the ob¬ 
servables lie within current experimental bounds. How¬ 
ever, upcoming experiments will be able to put stringent 
constraints on the TB and EB modes and a detection 
of one of these correlations could be the sign of a flavor- 
space locked gauge field. 

The perturbations in all three sectors are bound on all 
scales. We discover a secular instability in the scalar and 
vector perturbations of the theory, nonetheless, this is an 
artifact of a poor choice of initial conditions and does not 
describe a breakdown of the model. We analyze the effect 
of changing the YM coupling constant, the total amount 
of YM fluid and the initial distribution of energy in the 
electric and magnetic part of the fluid. Varying these 
initial conditions leads to big changes in the B cross and 
autocorrelation spectra. 

We compute constraints that future experiments put 
on the parameters of our theory assuming that TB and 
EB cross-correlation are measured. A detection of this 
coupling could be used to determine whether the gauge 
field is part of a dark sector that includes dark energy. 
We note that the effect of the flavor-space locked SU(2) 
electric or magnetic fields on particles with charge qym 
under the same SU(2) group would be negligible. How¬ 
ever, if dark energy couples to the rolling gauge field, or 
if the gauge field is dark energy, as in a gauge-flation sce¬ 
nario, then the rate of cosmic acceleration may be linked 
to the chiral asymmetry. 
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Appendix A: CAMB implementation 


Here we provide useful notes for implementing the 
Yang-Mills perturbations in GAME [30]. We will first 
outline the general approach to implement this model in 
CAMB and then mention some details in the subsections. 
Our universal approach is to compute the evolution of the 
background field ^(t) and then the perturbation equa¬ 
tions on top of these. The background field follows a 
simple ordinary differential equation (ODE) (Eqn. (8)), 
whereas the perturbations technically have to be solved 
as a system of PDFs due to their spatial variation. 

We solve Eqn. (8) with a fourth-order-Runge-Kutta 
method and store the results in an array. This happens 
in the subroutine init.background before CAMB ex¬ 
ecutes the actual Boltzmann integration. To put this 
model into the code we need to use the right units. 
CAMB is written in inverse Mpc units, and all energy 
densities are multiplied by 87rG. Therefore we define the 
background field in these units in order to not have to 
convert it to CAMB units everywhere, (j)Ci ^^nd the sub¬ 
script stands for CAMB. We dehne 


= 


Mp 


4>c = 


^Mp 


(Al) 


where (pQ and pc have mass units of 1 and 1/2 respec¬ 
tively. Physically this is a unusual redefinition, but nu¬ 
merically it makes the equations well behaved as we treat 
p and p' as two distinct quantities. With this mapping 
there is only one scale involved in the equations. Mp 
always refers to the reduced Planck mass Mp^ = SttG. 

There is a natural choice for initial conditions for the 
background held p. The energy density (Eqn. (7)) splits 
into the ‘electric’ (oc p') and ‘magnetic’ (oc gvMp'^) parts. 
Therefore, we introduce an angle 0 that captures the dis¬ 
tribution of energy between these two parts. Energy den¬ 
sity considerations can be used as initial conditions for p 
and p' by simply assuming a hxed fraction of the YM en¬ 
ergy density to the relativistic energy density. Inverting 
the relationship in Eqn. (7) gives the initial conditions 
(equivalent to Eqn. 11), 


P'c,i = 2 AYMf^7,o sin^ 9 


4>c,i — 


Ho 

gvM 


COS^ 6 


(A2) 


where we choose gyu, Rym and 0 G [0,7r/2], which con¬ 
trols the initial relative contributions of p and p', the 


electric and magnetic contributions to the YM energy 
density. Here H.y.o is the radiation energy density frac¬ 
tion today. For our analysis we choose 9 = 7r/4 to get 
a perturbation history close to vanilla ACDM. Deviating 
from this approach gives more exotic models that we ex¬ 
plored around Fig. 8. These initial conditions are called 
once before solving the background equation of motion. 

As shown in Sec. HI to first-order the YM fluid be¬ 
haves like radiation. For the energy density computa¬ 
tions in equations. f90 we can simply add a routine 
that multiplies the total radiation density by a factor of 
(1 -I- Rym)- The computation of the dark energy density 
has to be adjusted accordingly for the variable P7oOmegav 
in inidriver.f90. 

The Fourier transforms of the perturbation quantities 
are k dependent and should most efficiently be solved 
together with the vanilla ACDM perturbations in the re¬ 
spective routines: The ODE solver for CAMB is dverk 
which solves first-order ODEs. The second-order pertur¬ 
bation equations therefore have to be split into two first- 
order equations. The values of the functions and their 
derivatives get passed on to dverk by the subroutines 
derivs for the scalar and derivst for the tensor pertur¬ 
bations. This is where we modify the existing propaga¬ 
tion equations for the gravitational waves and add new 
expressions for the gauge field tensor and scalar pertur¬ 
bations. In the following sections, we will describe in 
some detail the implementation of this model in CAMB. 
The tensor modifications are more involved and have a 
bigger impact on the physics, and hence we will start 
with them. 


1. Tensors 

This subsection describes the implementation of ten¬ 
sor perturbation equations in CAMB. We will outline 
the general approach computing the Cps via the Boltz¬ 
mann integration and present the corresponding steps in 
CAMB. 

To obtain the temperature or polarization anisotropy 
for a given angular mode, one has to integrate the appro¬ 
priate transfer functions against the initial power pertur¬ 
bations from radiation domination onward until today. 
See Eqn. (9) of [53], using the appropriate initial condi¬ 
tions. The transfer functions are obtained by integrating 
the appropriate (tensor) source function against the ten¬ 
sor spherical eigenfunctions (Eqn. (15) of [53]). We need 
to alter the way the source function gets computed, as 
this contains all the underlying physics. 

Generally, the tensor perturbation quantities that 
enter the source function get defined in subroutine 
derivst and propagated in time in subroutine 
outputt. The routine derivst takes first-order ODEs 
and passes them on to the ODE solver used in CAMB, 
dverk. We add the new physics to derivst and pass the 
results on to the higher level ones. 

We extend the number of equations getting propa- 
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gated by 6 , 2 for the left handed gravitational waves 
(note, we turn one second-order ODE into two first- 
order ODEs), and 4 for the right- and left handed 
gauge field tensor perturbations tu/R- The variables 
getting passed on to dverk are labelled yt(l), yt(2), 

..., yt(n) and ytprime(l), ..., ytprime(n) for 
their derivatives. The corresponding physical quan¬ 
tities get mapped as described in Eqns. (A3). Here, 
n is the variable that has to be increased by 6 in 
subroutine SetupTensorArrayIndices. We are using 


the following prescription: in 

dervist and outputt: 

1 

II 

1 

y 

1 

II 


h'p = —k shear = —kY^ 

1 

1 

II 


hp = Hchi = Y 2 



1 

II 

y 

i;:_5 = -kY^-A 


II 

T" — hV' 

—^^n-2 

(A3) 

II 

1 

II 


II 

1 

Tr = F„_3 


1 

y 

II 

kY„_2 = y_3 


where we define 

Tl/r = i/2 a 

th/R 



to simplify the tensor propagation equation. Also, 
refers to ytprime (n) etc. 

In the following we will translate the perturbation 
equations into the language of CAMB. A dictionary for 
the tensor perturbations is given in Table II where we 
collect the most important quantities and procedures to 
identify the physics in the code more easily. We focus on 
the files cmbmain.f90 and equations.f90 which con¬ 
tains most of the physics involved in the Boltzmann in¬ 
tegration. Using the CAMB expressions, the left handed 
tensor perturbation equation (14) for becomes 


h'li -l- 2 — h'ji + + ^2^2 ~ ^ 


2V2 2V2 


2^2 


a 2 M|, 


'T'd - o I ,0 -I- „ gYM'p'^kTR 






= 0 


Y' 2 = —2 — shear 
” ^ a 


+ Hchi ^fc -I- {^YU^(i>c ~ 


H- ^(Pc^ri-2 ~ 2 


2Mp 


W-3 


/ gYM^ 


4- 


^/Mp 




where the operations written in red are changing sign 
going from right- to left handed propagation. Next, the 
gauge field tensor perturbation t in of Eqn. (14) becomes, 
using the background propagation equation ( 8 ) to elimi¬ 


nate the (t>'Q term 

Yn-2 =Yn-3 (^4-2gYM — k'j 

'Pc+9ym4'c'^ 

— V 2 (j)'Q shear 

where again the red operations change sign upon switch¬ 
ing from left- to right-hand polarized gravitational waves. 

For the perturbations, the initial conditions for the 
integration of the transfer functions against the initial 
power spectrum are defined in subroutine initialt. In 
vanilla CAMB only one tensor perturbation gets propa¬ 
gated. Hence we copy the initial conditions for the second 
gravitational wave history. 


- V2Hchi 



2. Scalars 

Here, we will collect results about the scalar imple¬ 
mentation in CAMB. See Table H for the definitions of 
the scalar variables in the code. The general strategy we 
will follow is analogous to the tensors: solving the two 
equations of motion for {Sm, S(j)} simultaneously. We will 
use the two constraint equations (Eqns. (17)) to solve for 
{y, y'}. The scalar perturbation equations that are going 
to source the anisotropic stress (Eqn. (18)) have to be 
translated to the synchronous gauge [43] as do the met¬ 
ric perturbation quantities b and $<3 from Eqns. (13) to 
be used inside CAMB. The expressions yield 


0 = a - 

n 

w v' 

<^G = V-^V-^-2na 




n 


where a = 2 ^ {k! -f In the calculations we need 
the first derivative of $< 3 . Computing this is intricate 
because it requires us to compute a'" and rj" which is 
not directly included in CAMB. However, this is actually 
not exotic: a description in a'" can be translated into one 
with the derivative of the equation-of-state parameter w' 
which is standard in dynamical dark energy models. The 
T derivative of 4 )g yields, in CAMB variables. 


^G = - 


aa'''ri 2aa''‘‘r] 


^/3 


a r] 
a' 


2sigmaa'' 


h''a 2sigmaa'' 


6 a' 


-b 


ka 

2sigmadota' 


ka"^ 


/csigmadota fcsigma 
3a' 3 


ka 

kz 

' y 


For 77 " we need h" and the CAMB variable sigmadot (see 
Table H). We use Eqn. (21c) from [43], which gives 


= -ij^SP-2m' + 2k‘^y 


(A5) 
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Theory 

Code 

cmbmain.f90 

kMkP,{k)\Aflf/ikr 

(Eqn. (9) of [53]) CalcScalCls 

Ahfk) computation 

y,Delta_p_l_k 

YtYi = d.TS^l^l{k,T)xi{ro - r) 

(Eqn. (15) in [53]) DoFlatIntegration 

equations 

.f90 ouput(t) and derivs(t) subroutines variables 

Source computation 

output(t) 

Eqn. (29) & (30) of [54]: -K.y & E 2 

Gaugeinterf aceEvolveTens: y(EVy,g_ix+2) & y (EVy,E_ix+2) 

Eqn. (29) of [54]: tTj 

outputt: pig for tight coupling 


opacity 

in [55] 

polter 

4'’ [55] 

pig 

Source tt in h" + 2^h' -|- k^h = n 

rhopi 

d/k from Eqn. (A6) and [43] 

vq 

e'/k 

ayprime (EVy,w_ix+l) 

h'/2k 

z 

density contrast 5 

clxq 

S' 

ayprime (EVy,w_ix) 

yk metric perturbation from [43] 

etak 

c? 

cs2 

TT^ = ^kTc{vb + a) (Eqn. (40) in [54]) 

pig = 32._dl/45/opacity*k*(sigma+vb) 

= Mp'^a^Y.iipi + Pi)vi 

= 2ky' dgq: total heat flux 

Mp'^a^ Yfi PiSi 

Total matter perturbation dgrho 

Yf,i Pi(^i 

Total Ma & Bertschinger [43] a dgs 

ky' 

ayprime(2) 

Y {h' + 6y') 

sigma 

i,ih'' + 6y”) 

sigmadot = -2*adotoa*sigma-dgs/k+etak 

v' 

1(sigma - z) 


TABLE II. Scalar perturbations: translations from theory to code. 


where we only need to plug in the expression for SP that 
is internally solved for in derivs 


SP = clSp+^ 


3— (1 + w) (cg — te) + w' 


(A6) 


where in our model we always have SP = Sp/3. Plugging 
this into the above Eqn. A5 gives the expression 

h” = -^5p-2Hh' + 2k'^p 

Ml 

which can be expressed purely in terms of GAME vari¬ 
ables. The expression for o'" can be obtained by taking 
a derivative of the second Friedmann equation and plug¬ 
ging the first one back in together with energy conserva¬ 
tion. The result reads 


a 


Iff 


( P w{3w-l) 4 

\Ml) 2V3 “ 


These conversions enter the scalar perturbation equa¬ 
tions. Also, we will have to convert the background field 
(p into the CAME variables as before. 


The equations need to be put into first- 
order form. We will use the prescription 

in the following equations to achieve this: 


Sp" = kYl 5m” = kYl_^ 

5(f)' = kYn 5m! = kYn -2 


5(fl = Yn-i 5m = Yn-3 

YLl = kYn K_3 = kYn-2. 

With these definitions the equations of motion for 5(f> 
and 5m become, where the index n refers to the extended 
set of equations in subroutine derivs. 


YI _2 = -kY„-i - kcp'cb + ky' 

Yi=l (29YMVc(i'n-3 - 3y„-,) - (.‘c (!=">' + i-t)) 


These equations get supplemented by the two constraint 
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equations that give y and y'\ 

^ fc(y„-2 - y„) + {2gW./M^b(t>l - 

y' = \/M^ ^(/)c {h' - $g) - 

+ 3(/>g& + 3V„-i — Yn-3 

which, as described above, get used to compute the scalar 
perturbation quantities 5m and 54>. 

This appendix serves as a guideline to implement mod¬ 
els like this one in CAMB without being too detail ori¬ 
ented. Following these descriptions results in relatively 
fast CAMB runs that take roughly four times the dura¬ 
tion of a vanilla run. 


Appendix B: Analytic solution to the background 
equation of motion 

The solution to the background gauge field equation of 
motion may be expressed in terms of Jacobi elliptic sine 
functions. Although these functions can be readily looked 
up in any math reference book, we give some of their 
properties here. To be succinct, the differential equation 

^+2f{tf=0 (Bl) 

has solution 

/(t) = cisn(cit-b C 2 I-1) (B2) 

where the function sn(u|m) is the Jacobi elliptic sine- 
amplitude function. This and related functions are de¬ 
scribed in Sec. 8.14 of Ref. [56]. The applicability of the 
above solutions to the description of both a scalar field 
(p with a Xcj)'^ self-interaction and SU(2) Yang-Mills was 
presented in Ref. [57]. 

A few useful facts are that sn(u|TO) can be evaluated 
by the following recipe: 

sn(u|m) = sine/), for u= . (B3) 

■JO yl — m sin^ 6* 

Our notation is consistent with Ref. [56] as well as Math- 
ematica [58]. (In some notation, m is replaced by 
in the integral on the right. Elsewhere, it is common to 
drop the m so that sn(u|m) becomes snit.) The deriva¬ 
tive is d/dusn(ulm) = cn(M|m)dn(u|m), which consists 
of Jacobi elliptic cosine-amplitude and delta-amplitude 
functions. Together they satisfy the identity 

sn(M|m)"^ -b cn(u|TO)^ dn(u|m)^ = 1. (B4) 


We are interested in the case m = —1, for which these 
functions oscillate with period T = r(l/4)^/-\/^. In the 
high-frequency limit, it may be useful to average over the 
oscillation period. In this case 


(sn(M| —I)^) = (dn(u| —I)^) — 1 = 1 — (cn(u| —1)^) 
87r2 

“ r(l/4)4’ 

(sn(u|-in = l. 

(B5) 

Illustrative plots of these functions are provided in 
Fig. 10. 



u 


FIG. 10. Jacobi elliptic functions sn(u| —1) and its derivative 
cn(u| —l)dn(u| —1) are illustrated. The dashed line in the 
first panel corresponds to a regular sine function with a scaled 
argument: sin(|u). 
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